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Abstract 

Extending recent numerical studies on two dimensional amorphous bodies, we characterize the 
approach of elastic continuum limit in three dimensional (weakly polydisperse) Lennard-Jones sys- 
tems. While performing a systematic finite-size analysis (for two different quench protocols) we in- 
vestigate the non-affme displacement field under external strain, the linear response to an external 
delta force and the low-frequency harmonic eigenmodes and their density distribution. Qualita- 
tively similar behavior is found as in two dimensions. We demonstrate that the classical elasticity 
description breaks down below an intermediate length scale £, which in our system is approximately 
23 molecular sizes. This length characterizes the correlations of the non-affine displacement field, 
the self- averaging of external noise with distance from the source and gives the lower wave length 
bound for the applicability of the classical eigenfrequency calculations. We trace back the "Boson- 
peak" of the density of eigenfrequencies (obtained from the velocity auto-correlation function) to 
the inhomogeneities on wave lengths smaller than £. 

PACS numbers: 46.25.-y Static elasticity, 72.80.Ng Disordered solids, 83.70.Fn Granular solids 
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I. INTRODUCTION. 

In a recent series of papers Q, 0, [J, we investigated the elastic response of zero tem- 
perature two dimensional (2D) amorphous systems. Our studies were motivated by the idea 
that such systems, although they appear perfectly homogeneous when looking at the density 
field, may be described as heterogeneous from the point of view of the theory of elasticity. 
The basic reason for this failure is now well identified: the underlying hypothesis of affinity of 
elastic deformations, implicit in standard elastic theory, needs not to apply to a disordered 
system. The relevant issue is therefore the scale above which a disordered, glassy system can 
be considered as homogeneous from an elastic point of view. 

Obviously, this question is important for the vibrational spectrum of such disordered 
systems; the excess of vibrational states at intermediate frequencies in the spectrum (the so 
called "boson peak" ) has previously been assigned to the existence of elastic heterogeneities 
J, ^]. Moreover, the field of plastic deformation of glassy materials, which has attracted 
considerable attention recently 0, Q, 0, ^, [ill [ill may be expected to be related to elastic 
heterogeneities. Other points of interest include the experimental evidence for dynamical 
heterogeneities in deeply supercooled systems [12|, which again could be expected to give 
rise to "frozen in" heterogeneities in low temperature systems. 

Our previous studies were limited to 2D systems, as this reduced dimension allows to carry 
out calculations on systems with large linear sizes using a limited number of particles. These 
studies allowed us to establish, for a standard computational model system, the existence of a 
length scale that can reach a few tens of particles, and below which classical elasticity breaks 
down. Similar conclusions were reached by Goldenberg and Goldhirsch Q]. This breakdown 
is revealed by a number of different diagnostics. Firstly, the so called Born expression for 
elastic constants is found to give incorrect results at zero temperature (where the fluctuation 
term does not contribute). This failure can be traced back to the importance of a non-affine 
contribution to the microscopic displacement field, while the derivation of the Born formula 
assumes affine displacement at all scales. The analysis of the correlation function of the 
non-affine contribution to the displacement field reveals a length scale £ over which this field 
is correlated, defining "soft" regions with large non-affine displacements. Second, the study 
of low frequency vibrations in these model disordered systems shows that the predictions 
of classical elastic theory are recovered only for wave lengths larger than £, meaning the 
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system is not homogeneous from the point of view of elastic properties below this scale. 
More recently, it was shown that the response to a point force is dominated by fluctuations 
for distances to the source smaller than £ j3j. Hence, £ characterises the self- averaging of 
the noisy response within each configuration which let us to call it the self- averaging length. 
The pressure dependence of £ has also been investigated in 2D demonstrating that the self- 
averaging length remains "mesoscopic" for low and moderate pressures, typically of the order 
of 40 particle sizes, but decreases at large pressures jl^j]. 

An obvious question that arises is the extent to which these results may depend, quali- 
tatively or quantitatively, on the dimensionality of space. Three dimensional (3D) systems, 
however, are considerably more difficult to study than the 2D case. In particular, if the 
elastic inhomogeneities take place on a length scale comparable to what is observed in 2D, 
very large systems have to be studied in order to reach the limit of elastically homogeneous 
systems, i.e. L ^> £ where L is the lateral size of a cubic system and £ the size of inho- 
mogeneities. Typical numbers of the order of at least 10 5 particles should be considered, if 
one wants to use the same tools and diagnostics in 2D and 3D systems. Although a number 
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of studies have appeared recently [la [la, [12J pointing to the existence of elastic inhomo- 



geneities in various types of disordered systems, all of them were realised for relatively small 
system sizes, making a direct comparison with our previous results difficult. In the same 
way, previous calculations of vibration modes in 3D systems have been limited to rather 
small sizes To our knowledge, this work is the first one to explore systems with 

lateral sizes that are appreciably larger than the expected scale of elastic heterogeneities. 

The aim of this work is therefore to characterise the elastic behavior of large 3D systems, 
using the same computer model and similar quench protocols as in our previous 2D studies. 
They are summarised in section [H] In section IIII1 we start by analysing the non-affine local 
displacement field in cubic samples submitted to uniaxial elastic deformation. From previous 
experience , we know that this type of analysis is the most cost effective in revealing the 
existence of inhomogeneities and their length scale. We then discuss the elastic response to 
a point force (section IIV[1 and corroborate why £ has been termed "self-averaging length" . 
Vibrational properties at very low eigenfrequencies - obtained by diagonalization of the 
dynamical matrix - are considered in section [V] and the density of eigenstates - computed 
by means of the finite temperature velocity auto-correlation function - in section IVII Our 
results are summarised in the last section. 
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II. DESCRIPTION OF SYSTEMS AND SIMULATION PROCEDURES. 



The initial configurations and their preparation are very similar to those described in 
Ref. fl for the 2D case. In order to make the comparison as direct as possible, the same 
type of potential was chosen, i.e. a slightly polydisperse Lennard- Jones potential Uij(r) = 
4e ((o"jj/r) 12 — (o"jj/r) 6 ) where the are taken uniformly distributed between 0.8<r and 1.2a, 
corresponding to a polydispersity index of 0.12. This is expected to be enough to destabilise 
a polydisperse crystal [21], and indeed no sign of crystallisation or demixing was observed 
in our simulations. The interaction energy scale e and the particle masses m will be taken 
to be strictly monodisperse. In the following, we will adopt the units appropriate for this 
Lennard- Jones systems, i.e. the mean diameter a will be our unit of length (and generically 
described as the "particle size"), while time will be expressed in units of r = a/ ma 2 /e. 

We studied various systems at constant density, p = N/L 3 = 0.98, which for our systems 
corresponds to a very low hydrostatic pressure at zero temperature (|P| ~ 0.2). The lateral 
size L of the system was varied between L = 8 and L = 64 (corresponding to N = 500 and 
N = 256000 particles). 

Disordered configurations are prepared by melting at high temperature (fcgT = 2e) 
an initially FCC configuration during 10 5 Molecular Dynamics Steps (MDS) using con- 
stant temperature molecular dynamics and velocity- Verlet integrator with a step size of 
O.OOlr. After the system was equilibrated, we checked that no crystalline order remains, 
and we begin the production run using two types of minimisation. The first one, called 
the "fast" quench, uses a direct conjugate gradient minimisation until (according to nu- 
merical tolerance) the zero temperature equilibrium state is reached. This type of minimi- 
sation was implemented for all system sizes. The second one, which we called the "slow" 
quench, is implemented for systems containing up to N = 62500 particles, corresponding 
to a lateral size L = 40. The protocol used in this case consists in equilibrating the liq- 
uid configuration at intermediate temperature (ksT = le) and then cooling this by stages 
{k B T = 5.10- 1 e, 10 _1 e, 5.1(T 2 e, 10" 2 e, 5.10" 3 e, 10" 3 e). At each stage, the system is "aged" 
(rather than equilibrated) during 10 5 MDS. Finally, the zero temperature state is reached 
using conjugate gradient minimisation from the last stage. Unless indicated otherwise, all 
the results discussed below refer to the fast quench procedure, which was carried out for all 
system sizes. 
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In order to obtain good quality statistics, this procedure was repeated between 1 and 10 
times, depending on the system size. The results presented in the following are averages over 
these different realisations of our amorphous systems. 
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FIG. 1: Non-affine part of the linear and reversible displacement field u(r) for the imposed macro- 
scopic uniaxial strain e xx = 10 -7 for a system containing N = 62500 beads (L = 40): (a) projection 
on the (x — z)-plane for all particles close to the plane. The length of the arrow is proportional to 
the displacement, (b) all beads of the same configuration with the 10% strongest non-affine dis- 
placements. (The short lines indicate beads with direct mutual interactions.) This subset of beads 
is strongly spatially correlated on short distances, however, it is homogeneously distributed and 
isotropic on larger scales. 

III. RESPONSE TO A MACROSCOPIC UNIAXIAL DEFORMATION 
A. Computational procedure and non-affine displacement fields 

In this section, we investigate the elastic behavior of zero temperature cubic samples, 
prepared as described above, submitted to an uniaxial traction. The procedure adopted 
is the following. First, a global deformation of strain e xx < 1 is imposed on the sample 
by rescaling all the coordinates in an affine manner. Starting from this affinely deformed 
configuration, the system is relaxed to the nearest energy minimum, keeping the shape of the 
simulation box constant. As a result, a displacement of the particles relative to the affinely 
deformed state is observed. This defines the non-affine displacement field u(r). 

A typical example for this field in the linear elastic response limit for a strain of e xx ~ 
10~ 7 ) is presented in the first panel of Fig. ^ It displays a two dimensional projection of 
u(r) in a plane containing the elongation direction for a system of size L = 40. (Note that 
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FIG. 2: Different moments of non-affine displacement field (y 2n ) 1 ^ 2n for n = 1, 2, 3,4 as a function 
of the imposed strain e xa: for systems of L = 40 obtained by means of the fast (open triangles) and 
the slow (full triangles) quench protocol. Both protocols show very similar results. The bold line 
on the left indicates the linear slope (<y 2n ) 1 / 2ra oc e xx . The vertical dashed line marks the limit of 
elastic response e p (L) ~ 10 -6 for L = 40. Also given is the residual plastic displacement field 
(obtained by reverse deformation back to the original macroscopic shape) for L = 40 and L = 32 
(full symbols). Residual fields below 10 -9 are due to numerical inaccuracies and the field can be 
considered as reversible. The sudden rise at e p for L = 40 corresponds nicely to the jump of the 
moments (y? n ) 1 ^ 2n . Note that the plasticity threshold e p depends strongly on the system size. 

projections on different planes are similar.) Visual inspection of such snapshots suggests 
that non-affine fields are strongly correlated over short and intermediate distances. This 
impression is also confirmed by Fig. ^b) where we focus on the 10% most mobile particles 
suggesting a connected cluster of these strong displacements spanning the simulation cell. 

In Fig. 121 we present the first moments ({u{e xx )) 2n ) l / 2n of the non-affine displacement 
field of system of size L = 40 averaged over all particles of an ensemble. Both quench 
protocols are very similar. (Only one moment is given for the slow protocol for clarity.) For 
e<£ p Ri 10~ 6 (vertical line) all moments are (up to prefactors of order one) identical which 
demonstrates an unique strain dependence for all beads. As one may also expect, a linear 
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strain dependence is found (bold line). At e p the moments increase suddenly and differ over 
more than an order of magnitude. This suggests an inhomogeneous strain dependence of 
the non-affine displacement field as will be discussed below (Fig. EJ). We note finally that 
the threshold e p decreases strongly with the system size (not shown). Hence, linear response 
requires much smaller deformations for large L. 

B. Plastic displacements and participation ratio 

The elastic (reversible) character of the deformations was checked by carrying out the 
reverse transformation and measuring the residual displacement of the particles, v_ i: which 
corresponds to plastic deformation. The moment of the residual field is indicated 

in Fig. El For e xx <C e p it is negligible and the deformation is, hence, elastic. Interestingly, 
elastic and linear elastic regimes coincide essentially as can be seen from the figure. For 
larger strains the residual displacements increase sharply over several orders of magnitude 
and coincide roughly with the n = 1 moment of the non-affine displacements. This shows 
that, for e xx > e p , the non-affine displacements are mainly due to plastic rearrangements. 

In view of the potential relationship with plastic deformation, it is interesting to inves- 
tigate in some detail the spatial features of the non-affine displacement field. Qualitatively, 
this can be achieved by representing, as shown in Fig. ^b), the particles that have the 10% 
largest non-affine displacements. This picture shows that the non-affine field is rather delo- 
calized, with the cluster formed by the most mobile particles spanning the entire simulation 
cell. A more quantitative view can be obtained by calculating the participation ratio for the 
non-affine displacement, defined by 



This participation ratio is shown in Fig. 01 for both quench protocols and L = 40. (A similar 
participation ratio may be calculated for the residual plastic field. However in this last case, 
the ratio at small e is due to numerics and at high strain it is identical to the participation 
ratio of the non-affine field.) Obviously, for sufficiently small deformations the displacements 
must depend identically for all beads on the applied strain and P has to become constant. 
As anticipated in Fig. El the presented data shows that this coincides with the linear elastic 
regime where all moments of the non-affine field are similar and the residual plastic field 
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FIG. 3: Participation ratio of the non-affine displacement field in a 3D system containing 62500 
particles (L = 40), as a function of the uniaxial strain e xx . Results for both fast and slow quench 
protocols averaged over eight and five configurations respectively. The given lines are guides to the 
eye. In the case of the fast protocol, a single configuration is also shown, for comparison with the 
averaged case. 



can be neglected. 1 The central point is here that the plateau value of the participation 
ratio is large (about 25%) indicating that the elastic non-affine displacement involves a 
substantial fraction of the particles. When the deformation exceeds the plastic threshold e p , 
however, the participation ratio falls rapidly, indicating that a plastic deformation proceeds 
via well localized events. 2 This is the first main result of this work. The implication from this 



1 Note that the jumps in the plateau value at small strains (below 10 -6 ) in Fig.0 specifically for the slow 
quench protocol, may be attributed to numerical inaccuracies due to the small displacement fields which 
have been compared. 

2 It should be stressed that the more gradual decay of the participation ratio, in the fast protocol case 
(Fig[3J), is due to ensemble average. It describes essentially the probability that no jump has occurred for 
smaller e values, and it is thus related to the distribution of the plastic thresholds e p from one configuration 
to the other. In the case of a single configuration, the decay of the participation ratio is more sudden, 
and can be compared with the averaged decay in the slow protocol case, where the distribution of e p 
is narrower. This is, in fact, the main difference observed between the two quench protocols: the plastic 
threshold e p is well defined in the slow quench protocol case, but it is largely distributed (over more that 
one decade) in the fast quench protocol case. Note however that we do not have, at this stage, a sufficient 
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difference in behavior is that the localized events occurring in plastic deformation cannot 
be directly inferred from the general pattern of non-affine displacements. This does not 
mean that plastic displacements and strong non-affine elastic displacements are completely 
uncorrected. 

In other words, energy barriers (which are relevant for plastic deformation) are not directly 



24j . Interestingly, the main influence of 



related to the local curvature of the energy minima 
performing a slow quench seems to be that the plasticity limit is increased, meaning that the 
system has been brought to a slightly more stable configuration with higher energy barriers 
without, apparently, changing measurably the local curvature of the energy minima. In fact, 
properties such as the vibrational modes discussed below, are much less affected by the 
quench protocol. 

In the reminder of this paper, we focus on the linear elastic response. We normalise the 
non-affine field by its second moment, i.e., is replaced by uj (u 2 ) 1 ^ 2 , in order to consider 
a strain independent reduced field. 

C. Hydrodynamic limit: Lame coefficients 

We turn now to the calculation of the Lame coefficients A and fi, which characterize 
the elastic behavior of an isotropic medium in 3D (22 1. Our results for these coefficients as 
a function of system size are shown in Fig. 0] This figure compares two different ways of 
obtaining the coefficients. A a and fi a are obtained under the assumption that the non-affine 
contribution to the total displacement field of the particles is negligible. They are simply 
the Born estimates, that can be, in a system with pairwise interactions, computed from 
the reference configuration by carrying out a simple summation over all pairs of interacting 



particles (see for example ref. p) 

A a = * = ± E (W*) - ^'M) ^ ( 2 ) 



where U is the interaction potential. The second estimate corresponds to the "true" value 
of the elastic coefficients, obtained by computing the actual stress in the sample after the 
relaxation that introduces the non-affine part of the displacement field. The averaged stress 



number of configurations to compute more precisely this distribution. 
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FIG. 4: Lame coefficients A (spheres) and fj, (squares) vs. system size L. Full symbols correspond to 
the direct measurement using Hooke's law, open symbols are obtained from Eq. (j2J which supposes 
affine deformations (Born term). The effect of system size is weak. For large boxes we get /i ~ 15 
and A ~ 47. The coefficients relying on a negligible non-affine field differ by a factor as large as 
2 from the true ones. Clearly, a calculation taking into account the non-affine character of the 
displacement is necessary for disordered systems. 

is computed using the usual microscopic expression 

°<*f> = J3 ^2(Zij)a(Fij)f> (3) 

i,3 

where the Greek indices refer to cartesian coordinates, and is the force between particles 
% and j. 

The Lame coefficients are then obtained from the standard formulae a xx = (X+2/j 1 )e xx and 
a yy = Xe xx for a deformation tensor which has only an e xx component (e xx is here the global 
deformation imposed on the sample). For larger systems, we obtain /j « 15 and A ~ 47. 
Hence, we find that the true values of A and \x differ considerably from the Born estimates 
which indicates the importance of non-affine displacements in determining the stresses in 
the material. This contribution tends to lower the shear modulus /x, and to increase the 
coefficient A. From the measured values of A and /i we get a bulk compression modulus 
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K = A + 2/i/d ~ 57, a Young modulus E = 8Kfi/(3K + /i) ^37 and a Poisson ratio 
v = (3K — 2{i)/2(3K + /i) « 0.4 22j. Remarkably, the bulk modulus would be correctly 
predicted by the Born calculation. This means that the non-affine part of the deformation 
does not contribute significantly to the increment in the isotropic pressure under compression 
or traction, but is mainly associated with shear deformations. (This point will be further 
elucidated below when we will discuss Fig. 0) Such a situation would be natural in high 
pressure systems, in which the repulsive inverse power part of the potential dominates the 
interaction, and compression can be accommodated by an affine rescaling of all coordinates. 
It is, however, less expected in our low pressure systems. 

D. Correlations in the non-affine displacement field 

The preceding results call for a more thorough analysis of the correlations of the non- 
affine displacement field which apparently can not be neglected for macroscopic quantities 
and should therefore be even more relevant for finite wave length properties. 

Following Ref. , the non-affine correlation field can be analyzed by computing the 
correlation function C(r) = (u(r) ■ u(0)). (The averages are taken over all pairs of monomers 
(i, j) being a distance r apart.) As can be seen in Fig. a decay over a typical length of 
23 particle sizes is observed (bold line), before the correlation function exhibits a negative 
tail (see inset). The 2D case (disks) is also included for comparison. Qualitatively similar 
behavior is found. The anti-correlation can be associated visually in 2D with the solenoidal 
character of the non-affine field . The organisation of the non-affine deformation in 
"vortices" is less obvious in 3D (see inset) as manifested by the about seven times weaker 
amplitude of the negative tail. However, this description can be compared with the direct 
visualization of Fig. ^ where the intricate structure of the vortices is shown. See also Fig. 
below. 

That the displacement field is indeed correlated over a size £ ^> a is further elucidated 
in Fig. |H1 Here we consider the systematic coarse-graining of the displacement field 



^#) = ]v-5>fc) ( 4 ) 

of all Nj beads contained within the cubic volume element Vj of linear size b. In the figure, 
we have plotted the (normalized) correlation function B(b) = ((U_j(b) 2 ) j / (u 2 )) 1 ^ 2 versus the 
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FIG. 5: Correlation function C{r) of the non-affine deformation field as a function of the distance 
between pairs of beads r. The correlation function is averaged over all possible pairs. The full line 
corresponds to a 3D sample containing 256000 particles (L = 64). Note that the "fast" quench 
protocol and the "slow" quench protocol (not shown here) give the same curve. Data from for 
a 2D system of linear size L = 104 is shown for comparison (line with open dots). In the 2D case, 
results are averaged over 30 configurations. The data are plotted in linear scale. Inset: Enlargement 
of the anti-correlation regime. Note the negative - although weak - correlation of the 3D correlation 
function for r > 23a. 



size of the coarse-graining volume element b (normalized by L). Two and three dimensional 
systems are again compared. In both cases we find an exponential decay which is well fitted 
by the characteristic scales £ ~ 23 for 3D and £ « 42 for 2D. Apparently, £ is similar to the 
distance where C(r) becomes anti-correlated. Note that the total non-affme displacement 
field of the box must vanish - since the center of mass of the system is fixed- and therefore 
B(b) — > for b/L — ► 1. This sum rule explains the curvature in the data and the sharp 
cut-off on the right hand side of the figure. The nice agreement of both estimations of £ 
demonstrated in the Figs. 03 and El and the similar behavior in both dimensions is the second 
central result of this work. 

More systematically, the displacement field can be investigated in Fourier space U_(k) = 
YliLi n(ju) ex P(ik ' Li)- Apart the normalisation factor 1/Nj, this is close to the volume 
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FIG. 6: The (normalized) magnitude B(b) of the non-affine field averaged over a volume element 
of lateral size b is traced as a function of b/L. The plot confirms that the correlations decay with a 
characteristic length £ = 23 for 3D and £ = 42 for 2D (spheres, for 10 configurations) respectively. 

element coarse-graining method, Eq. (J1J). As usual, the fluctuations can be directly computed 
from S(k) = ^p(||L[(^i)|| 2 ) where the average is taken over all wave vectors with k = \ \k\\. Our 
results for S(k) are not presented since we prefer to discuss below in Fig. the longitudinal 
and transverse contributions, Si{k) and Sx(k). Note also, that S(k) is related to the real 
space coarse-graining correlation function B(b) of linear size b = 27r/k given in Fig. |3 3 
Interestingly, it can be readily shown that 

S(k) = — exp(ik ■ (r { — rj)(u(r i )u(r : j)) = (u 2 ) + f drexp(ik-r)C(r)pg(r) (5) 



3 The wave vector k must be chosen commensurate with the periodic simulation box. For this reason fewer 
values of the wave length A = 27r/||fc|| = L/ \f{n 2 +m 2 + l 2 ) (n, m and / being integers) can be computed. 
This is the principle disadvantage of the Fourier transformation of the non-affine field for computational 
studies in 3D - where the system size is necessarily rather limited - compared to the real space coarse- 
graining method which allows also a continuous variation of the box length b. Moreover, it is easy to show 
that B 2 (b) is the average over a volume b d of the correlation function C(r) whose Fourier transform is 
related to S(k). Thus a clearly shown size dependence of B(b) can be masked in the power law decay of 
S(k) (at large k). 
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FIG. 7: The squared amplitudes of the Fourier transforms k 2 Si and fc 2 Sr for the div (lower data) 
and curl (upper data) of the non-affine deformation field (see Eq. ([?))) plotted as functions of the 
wave length A = 27r/||fc||. Different system sizes, as indicated in the figure, have been included to 
demonstrate that Si and St are essentially system size independent for small A. Note that the 
statistics deteriorates for large A. 

with C(r) being the displacement correlation function discussed in Fig. and g(r) the stan- 
dard pair distribution function. Since g(r) becomes rapidly constant it follows that S(k) cor- 
responds essentially to the Fourier transform of C(r). This demonstrates that (very roughly) 
C(r) and B(b) (shown in Fig. EJ and |B1 respectively) are related by Fourier transformations 
and explains why similar characteristic sizes £ 3> u have been obtained. 

We demonstrate now that the non-affine displacement field in 3D is indeed of predomi- 
nantly solenoidal nature. This has been anticipated by our previous studies on 2D systems 
jj| and by the values of the elastic moduli fi and K discussed in Sec. IIII CI The transverse 
and longitudinal contributions to the displacement field can be numerically obtained by 
computing 

Z(k) = ~k A (kA U(k)) , L{h) = TiHk- Uik)) (6) 
k z k z 

Obviously, U = T + L, k-U = k-L, kAU = kAT and k-T = kAL = 0. The norms of 
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these quantities, for instance 

k 2 S T (k) = ^(\\T(k)\\ 2 ) = ^< \\J2^ A ^ e Mih-r U )\\ 2 > (7) 

i 

and similarly for the longitudinal part k 2 Si(k), are the Fourier transforms of V A u and 
Vu . They are plotted in Fig. [7| as functions of A = 2ir/k. We find that the longitudinal 
contribution (data points at the bottom) is about 10 times smaller than the transverse one. 
Note that all data points of different system sizes collapse well for wavelengths corresponding 
to the non-affine displacement regime, A < £. We find that k 2 ST{k) is more or less constant 
while k 2 SL(k) decreases weakly. Since T and L are orthogonal this yields S(k) = 5r(&) + 
Si,{k) ~ Sr(k) and, hence, the algebraic relation S(k) oc 1/k 2 for large k. 4 Unfortunately, 
for larger A the statistics deteriorates due to the smaller number of wave vectors which can 
be considered. Our data may suggest that Sr{k) increases for A > £, however, new data 
with larger boxes and with better statistics is warranted to confirm this. 



4 see previous footnote. 
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FIG. 8: Comparison of the incremental stress fluctuations bo a p = (< <r£g > — < a a p > 2 )*/ 2 with 
the mean vertical component of normal stress < a zz > for the ensemble average discussed in text, 
measured along the vertical line [x = y = 0). The mean stress decreases essentially as 1/V 2 , as 
expected (without logarithmic corrections) in 3D for positions far from the source and the fixed 
walls. (Note that < a zz > has to decrease less rapidly close to the walls, r ~ L/2.) In contrast, the 
fluctuations decay exponentially over the whole available system. The characteristic length scale is 
similar to the size £ ~ 23a of the correlated non-affine displacements. 

IV. SELF-AVERAGING OF THE RESPONSE TO A POINT FORCE 

In 2D, we showed previously j^| that the deviations from continuum elasticity at small 
scales could be revealed by studying the response of the system to a localised force. The 
same effect is illustrated for 3D in Fig. |H1 This plot is obtained as follows. A small force 
is applied to the particles contained in a small region of space (sphere of diameter 4). The 
force is applied in the z direction, and its magnitude is chosen small enough to remain in 
the (linear) elastic region (F z = 10). In order to maintain global force balance, the system 
has periodic boundary conditions in the x and y directions, but is immobilised by two fixed 
walls in the z direction. The dimensions of the simulation cell, which contains N = 165000 
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particles, are L z ~ 105 in the z direction, L x = L y ~ 40 in lateral directions. The stress 
tensor is then measured within the sample using the standard Kirkwood definition 
calculated in small rectangular boxes of fixed size (3, 3, 5) centered in (x, y, z). Those boxes 
are displaced in all the material by unit steps of one in the three directions. For each step the 
six components of the stress tensor are calculated and averaged on a statistical ensemble of 
200 configurations. Such an ensemble is obtained by taking 10 independent configurations, 
and for each configuration by changing the position of the point source within the sample 
reindexing the configuration in such a way that the origin of the source still is placed in the 
midplane, equally distant from the two fixed walls. In the absence of the external force, the 
local stresses in amorphous systems are usually nonzero ("quenched stresses"). The relevant 
quantity that defines the response to an external force are therefore the incremental, rather 
than total, local stress tensor. Once such a calculation done, the quantity of interest is the 
fluctuations of the stress tensor in the statistical ensemble j3j. 

In Fig. both the average value and the fluctuations of the stress tensor are shown. The 
average response is compared to the prediction of continuum elasticity, which appears to be 
perfectly obeyed on average even very close to the source. This average response exhibits 
the 1/r 2 decay characteristic of the Green function of classical elasticity in 3D. However, up 
to length scales of 50 the fluctuations are considerably larger than the average value of the 
stress. The fluctuations, on the other hand decay exponentially away from the source, with a 
characteristic length 23 - the same value as obtained from the correlation functions discussed 
in the previous section. In fact, the relative stress fluctuations, for instance d~cr zz / (o~ zz ), scale 
like exp(— r/£)r 2 and show non-monotoneous behaviour (not shown). They increase first up 
to 2£ (due to the decreasing mean stress), but decrease ultimately exponentially. With other 
words, self-averaging (within every configuration) of the noisy signal occurs on distances set 
by £ which we call, hence, the self- averaging length. 5 Unfortunately, our simulation boxes 
are yet too small to illustrate more clearly the exponential decay of the relative noise for 
r > 2£ » 50. 



5 Absolutely similar behaviour has also been found in 2D systems |2|. There the relative fluctuations scale 
like exp(— r/£)r, i.e. exponential screening sets in for distances r 3> £. Since £ is about twice larger in 
2D it turns out that about the same linear distance is required in both dimensions for self-averaging to 
become effective. This suggests that if different dimensions d are compared one should rather use the 
notion "self-averaging length" for (d — 1)£. 
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FIG. 9: Eigenfrequencies uj of the first 100 modes as a function of mode number p. The frequencies 
are rescaled using the system size and the sound velocity (obtained from the Lame coefficient 
in Fig. 0J). The horizontal lines correspond to the results expected from macroscopic elasticity. 
Frequencies of small systems are systematically too low. Besides the obvious Goldstone modes 
(j) = 1, 2, 3), all frequencies are finite. 

V. LOW FREQUENCY EQUILIBRIUM VIBRATION MODES 



We turn now to the determination of low frequency vibration modes. These were deter- 
mined from a direct dia gon alization of the dynamical matrix, using a modified version of 
the PARPACK package |23|. For a periodic cubic system described by classical elasticity, 
the structure of the low frequency end of the spectrum is well known. Each mode is char- 
acterised by a wave vector k = ^f(l,m,n). Transverse (resp. longitudinal) have frequencies 



UJ 



cxk (resp. u = C£&), where the sound velocity is given by ct = y/jijp ~ 4.2 (resp. 
cl = a/(A + 2/i)/p « 8.9). As a result, the modes should have well defined degeneracies. For 
example, the lowest lying mode (±1, 0, 0) should have 12 fold degeneracy, corresponding to 
the two transverse polarisations for the 6 wave vectors of length 2tt/L. The second frequency 
has degeneracy 24, and so on. In our previous analysis of 2D systems [l], we found that this 
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degeneracy of the low frequency modes was lifted for small systems sizes. A scaling analysis 
of the modes allowed to establish in a different manner the length scale £, above which the 
medium can be considered as elastically homogeneous. 

Our results for the low frequency modes of three dimensional systems are shown in Fig. El 
Plotting the rescaled and averaged frequencies (y) = ((ujL/2ttct) 2 ) as a function of mode 
number p, 6 it is clear that the degeneracies and the associated step-like behavior begin 
to show up only for large systems, containing at least 32000 particles (lateral size 32). For 
the largest systems however (lateral size 64) the discreteness of the low frequency spectrum 
is well apparent, typically up to the 4th eigenfrequencies. In view of the large value of 
cl compared to ct {cl — 2.1cr in our system) we have concentrated on the analysis of 
transverse modes. Longitudinal modes enter only at higher frequencies, and are mixed with 
shorter wave length transverse modes, making their contribution more difficult to identify. 
If we take as a criterion the existence of a gap separating the first 12 eigenfrequencies from 
the rest of the spectrum, it appears that the minimum size for applying continuum elasticity 
is comprised between L = 16 and L = 32 particle sizes. 

This analysis can be refined using a scaling plot of the mode frequencies as a function 
of the "theoretical" wave length, or more precisely the wave length of the elastic wave 
that would appear in the spectrum with this mode number according to elastic theory. 
Fig. El is constructed by averaging, for each size, the frequencies that correspond to the 
first elastic mode in elastic theory (e.g. the first 12 frequencies are averaged to obtain the 
lowest frequency point, the next 24 for the second point, and so on). The resulting frequency, 
divided by the value expected from elastic theory, is plotted as a function of wave length. 
Note that all data points collapse on the same master curve irrespective of the box size L. 
The plot shows that when the wave length is lower than the self-averaging length, deviations 
from elastic theory become significant, whatever the size of the system. This estimate for 
the size of elastic inhomogeneity is therefore in fair agreement with those obtained in two 
previous sections from the analysis of the linear response to an external load. 



6 The rescaling eliminates trivial dependencies on system size and sound velocities which are anyway ex- 
pected from continuum theory. The rescaled frequencies are averaged over the ensemble of available con- 
figurations. Note that differences between the eigenfrequencies of different configurations of the ensemble 
are weak, however, even for small system sizes. 
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FIG. 10: Eigenfrequencies corresponding to the first four levels predicted by continuum elasticity, 
plotted as a function of the corresponding wavelength A = 2tt/ || k ||. The eigenfrequencies uj have 
been rescaled by their continuum theory values. This yields a perfect collapse of all data points. See 
main text for details. The crossover to continuum behavior (indicated by the horizontal line) takes 
place, as expected, at the self- averaging length A ~ £. The frequencies decrease systematically with 
smaller A. Inset: The same plot for the slow quench protocol shows barely distinguishable behavior 
(same symbols as in the main plot). 
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FIG. 11: Density of states g(u>) for 3D amorphous systems of different system sizes L. The lines 
indicate three power law slopes, the dashed one being the Debye prediction g_o(u>) calculated 
from the known sound velocities. The dashed-dotted linear relation corresponds to the non-affine 
displacement field regime (A <C £). Also given are the characteristic frequencies ujl,t = cl,t^ '/£ 
associated with the self-averaging length £ and the Debye frequency u>d ~ 18.3. Larger frequencies 
correspond, in fact, to vibrations on very small scales, A <C o. Inset: g{uj)/uj 2 vs. lo for L = 56 
(same symbols as main figure). Note that ojt and col describe correctly the position and width of 
the boson peak. 

VI. DENSITY OF VIBRATIONAL STATES 

The (normalized) density of vibrational states (DOS) of a 3D solid may be denned by 
g(uj) = jrjy X/p=i fi( u ~ u p) w ith co p being the harmonic eigenfrequency corresponding to 
the mode number p. Hence, for small systems (of order of 10 3 beads) one can compute the 
complete DOS from the eigenfrequencies extracted by exact diagonalization of the dynamical 
matrix, just as we have done in the previous section. Obviously, for systems containing about 
10 5 particles the number of modes one may compute is rather limited. From the 100 modes 
we have presented in Fig. 01 one estimates roughly oc p a with a ~ 1. Hence, the DOS 
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increases approximately linearly, g(uA oc uj 2 l a ~ 1 ~ uj, for small uj. 

Following standard procedures we have instead obtained g{oS) by Fourier transfor- 
mation of the velocity auto-correlation (v(t) ■ v(0)). In contrast to the previous sections, we 
consider here configurations at finite, yet very low temperatures T. For the data presented 
in Fig. we have used T = which is three orders of magnitudes below the glass tran- 
sition of our systems. We start with quenched configurations at T = which we subject to 
a Maxwell velocity distribution. Following a thermalization phase of St = 10 3 the velocity 
correlation function is sampled over St = 100. Different temperatures have been checked 
and we have verified that the DOS becomes temperature independent at low T (not shown). 
As can be seen from the Fig. El our results become rapidly system size independent for 
sufficiently large frequencies uj ^> ct,l2vt/L. We remember that only uj values corresponding 
to wave vectors k compatible with the box size are physically acceptable in the continuum 
limit. The corresponding finite size effects at low frequencies can be seen on the left hand 
side of the figure. 

In linear coordinates, g{uj) is roughly symmetric around its maximum at uj ~ 14.3 and 
may be very crudely described as linear for small uj in agreement with the estimate from Fig. El 
described above and the dashed-dotted line indicated in the main panel of Fig. ^2 Note that 
the maximum is slightly smaller than the Debye frequency ujd = (18p7r 2 /(l/c 3 - / + 2/cfr)) ~ 

18.3 for our systems. (The Debye frequency is in turn smaller than the frequency ct27t/ct" « 

26.4 associated with a wave length of monomer size.) The log- log plot presented in the main 
figure shows various frequency regimes. For very small frequencies our data is roughly in 
agreement with the Debye continuum prediction goi^) — 3uj 2 /ujp (dashed line). The DOS 
increases than more rapidly with frequency up to ujt = ct2tt/^ m 1.1 - corresponding to a 
wave vector given by the self-averaging length - where g(uj) has power law slope of exponent 
2 (bold line). This can be more clearly seen in the inset featuring the enigmatic Boson peak. 
Apparently, the width of this peak is well described by u T and the frequency ul ~ 2.3 for 
the corresponding longitudinal wave. Hence, the boson-peak is fixed by the self-averaging 
length and marks the crossover between the continuum elastic behavior (dashed line) and the 
non-affine displacement field regime, where g(oj) oc uj (dashed-dotted line), at larger uj and 
smaller wave length A. Since continuum theory overestimates the frequencies for A £ (see 
also Fig. ITTH) this implies an excess of modes at smaller frequencies. Apparently, these modes 
are shifted to the edge of the non-affine regime. The boson-peak is merely a consequence of 
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the inapplicability of the continuum theory at A <C £ ~ 23. 



VII. CONCLUSION 

We have investigated the approach of the continuum limit for elastic properties of 3D 
amorphous systems and compared our computational results with our previous work on 
similar 2D systems. The results are extremely similar in both cases, and can be summarized 
as follows. The elastic constants estimated using the Born formulae are not accurate even 
at zero temperature, therefore revealing the importance of the non-affine component of the 
deformation field. This non-affine deformation field, which affects mostly the shear response 
(as compared to compressibility) is correlated over intermediate distances, of the order of 



ength is significantly smaller than 
17| . but implies that rather large 



23 interatomic distances in our case. This correlation 
in 2D, in agreement with the findings of Rossi et al 
samples should be used to discuss elastic or vibrational properties of 3D systems as well. 
By considering the Fourier transformation of solenoidal and longitudinal part of the non- 
affine field we have demonstrated (Fig. |7J) that the 3D non-affine field is mostly rotational 
in nature, in agreement with the visual impression of snapshots. The response to a delta 
force perturbation allowed us to measure the self-averaging of the noisy response within 
a configuration. The stress fluctuations decay exponentially with distance from the source 
with a characteristic length given by the correlation length £ of the non-affine field which 
suggests to us the notion "self-averaging length" . 

Vibrational modes are obviously strongly affected by the existence of elastic hetero- 
geneities, and cannot be predicted using elastic theory if their wave length is too small. 
From our scaling analysis (Fig. H"0|) it appears that the frequencies are smaller than ex- 
pected from continuum theory, therefore implying an excess of modes in the low frequency 
region compared to the Debye prediction. This excess has been analyzed in Fig. ^2 showing 
the density of vibrational states and demonstrating that the "boson peak" is located at the 
edge of the non-affine displacement field. That both position and width of the peak are given 
by the self-averaging length £ is the central novel result of this work which has not been 
previously established for 2D systems. 

The focus of this work has been primarily on the linear elastic behavior of amorphous 
solids. Our preliminary study of larger (uniaxial) deformations that go beyond the elastic 
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limit indicates that plastic events are rather localized individual events characterized by a 
very low participation ratio. In the recent work |15] de Pablo and coworkers pointed out 
the possibility of regions of negative shear modulus in quenched amorphous systems - such 
regions being stabilized by the "normal" material in which they are embedded. The typical 
size of these regions is much smaller than the size for elastic inhomogeneities discussed in 
this work, implying they are more likely to be linked to elementary rearrangements taking 
place at the onset of plastic deformation, which usually imply small numbers of particles 
even localization along a shear band [lOj, |Tjj . Such a difference in elastic 
and plastic deformations was also observed for 2D systems in Ref. 

The general picture that emerges is therefore that of a hierarchy of length scales. Disorder 
at the level of 2-3 atomic distances can be interpreted as implying the existence of regions 
with negative moduli, which will give rise to plastic yield. On a larger scale, this disorder gives 
rise to strong non-affine displacement fields in elementary deformations. Finally, convergence 
to standard continuum properties is obtained over length scales larger than the self-averaging 
length. In our analysis, carried out for a typical liquid state density and at zero-temperature, 
£ is found to be large, but finite. In analogy with what is found in 2D, we expect it to decrease 
with increasing density, and possibly to diver ge a s the density is lowered and the system 
loses mechanical stability, as suggested in Ref. |l6| |. 
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